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Abstract. In this paper, we study the effect of MHD turbulence on the dynamics of dust particles in protoplanetary 
disks. We vary the size of the particles and relate the dust evolution to the turbulent velocity fluctuations. We 
performed numerical simulations using two Eulerian MHD codes, both based on finite difference techniques: ZEUS- 
3D and NIRVANA. These were local shearing box simulations incorporating vertical stratification. Both ideal and 
non ideal MHD simulations with midplane dead zones were carried out. The codes were extended to incorporate 
different models for the dust as an additional fluid component. Good agreement between results obtained using 
the different approaches was obtained. The simulations show that a thin layer of very small dust particles is 
diffusively spread over the full vertical extent of the disk. We show that a simple description obtained using the 
diffusion equation with a diffusion coefficient simply expressed in terms of the velocity correlations accurately 
matches the results. Dust settling starts to become apparent for particle sizes of the order of 1 to 10 centimeters 
for which the gas begins to decouple in a standard solar nebula model at 5.2AU. However, for particles which 
are 10 centimeters in size, complete settling toward a very thin midplane layer is prevented by turbulent motions 
within the disk, even in the presence of a midplane dead zone of significant size. These results indicate that, when 
present, MHD turbulence affects dust dynamics in protoplanetary disks. We find that the evolution and settling 
of the dust can be accurately modelled using an advection diffusion equation that incorporates vertical settling. 
The value of the diffusion coefficient can be calculated from the turbulent velocity field when that is known for a 
time of several local orbits. 
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?H ■ 1. Introduction planet is to form subsequently. Drag is responsible for 

I example for the radial migration of solid bodies toward 
Dust particles are very likely to be the basic building -^^^^ ^^^^^ ^ ^-^^ j^^, ^^-^^ pressure decreases 
blocks that need to be assembled to make planets. ^^^^^^y outwards. The timescale for this migration is so 
In the "core-accretion" model for the formation of gi- ^apid for centimeter and meter size bodies that it has be- 
ant planets (|Mizuno|l980t [Pollack et aljl996|), the growth come a problem for s tandard planet formation theories 
of their solid cores proceeds through the accretion of ob- dWeidenschillind Il977l) . a resolution of which might in- 
jects ranging from micron sized dust particles to planetes- ^olve the particles accumulating close to pressure maxima 
imals of radius ^10 km to eventually reach about ISM^. possibly at the center of vortices 
Dust is also one of the main observational tracers of the 

struc t ure of protoplanetary accretion disks ijAdams et al.l Another important aspect of dust dynamics is the ten- 

Il987t iD'Alessio et allEool . A detailed knowledge of its dency to settle towards the midplane of the disk, which in- 

dynamics is therefore needed both in order to make the- creases with particle size. Because of this , centimeter and 

oretical models and to give the best interpretation of the metre sized bo dies accumulate clo s e to th e disk midplane. 

observations. As argued bv iGoldreich fc WardI l)l973|) . a gravitational 

The key ingredient that drives the evolution of instability in the dust sub-disk in triggered when the dust 

dust particles is the d rag force they feel from the gas density becomes high enough fse e alsoiSafronov 1969, and. 

l)Weidenschilling| Il977|l . Gas must be present if a giant as remarked in a note in proof bv lOoldreich fc Wardll973l 

iGiirevich fc Lebedinsk vll95(ll . At this point, larger bodies 

Send offprint requests to: S.Fromang known as planetesimals form. 
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However, because the pressure force acts only on the 
gas , a velocity shear between the dust sub-disk and the 
upper layers of the gas disk develops. This can lead to the 
so called "shear-induced" turbulence that may be able 
to mix the dust layer enough to prevent enough gravi- 
tational set thng to satisfy the condition for gravitationa l 
instability l)Garaud et alJl2004l: Icdmez fc OstrikeJl2005|) . 
However, the onset of this shear induced turbulence de- 
pends on the vertical profile of the dust sub-disk which is 
quite difficult to calculate. Indeed, this is determined as 
the result of the interplay between various physical pro- 
cesses, such as d ust coagulation, B rownian motions and 
turbulent mixing ijPuUemond fc Dominik..2005.) . The pur- 
pose of this paper is to investigate the behaviour of the 
dust disk when just one of these operates, namely MHD 
turbulence, from first principles. 

Turbulence in protoplanetary disks is believed to be 
the result of the magnetorotational instability (MRI). 
More than a decade of analytical an d numerical work (see 
a review bv lBalbus fc Hawlevlll998il has shown that this 
is capable of leading to sustained turbulence, and mag- 
netic field with associated stresses that transport angu- 
lar momentum outwards. More and more of its proper- 
ties have been analysed thanks to the increased compu- 
tational power currently available. However, little work 
has focus ed on the effects of MH D t urbulence on dust 
dynam ics. ICarballido et all l)2005j) and lJohansen fc Klahil 
l|2005|) analysed the radial diffusion of small particles 
and calculated the effective Schmidt number, which is 
the ratio between the anomalous turbulent viscosity and 
the anomalous diffusion coefficient. They found differ- 
ent result s which could be du e to their different ini- 
tial se tup. Ijohansen et al.1 1)2005|) and iFromang fc Nelsonl 
1)2005(1 studied the radial migration of centimeter and me- 
ter size bodies in turbulent disk models and found that 
dust trapping in local pressure maxi ma could be an impor- 
tant confining effect. Very recently, iTiirner et al.l l|200dl 
investigated the dispersion of a passive scalar in stratified 
disk models by releasing particles at high altitude in the 
disk. They interpreted the effect of MHD turbulence in 
terms of a damped wave equation. However their analy- 
sis was restricted to infinitely small particles. They did 
not consider the vertical settling that begins to occur for 
larger particles. Thus in their case stratification is only sig- 
nificant for the gas, being felt by the dust particles only 
through their strong coupling to the gas. 

A first attempt to study vertical settling numer- 
ically was made through g lobal SPH calculations by 
iBarriere-Fouchet et al.l l|2005|) . They confirmed earlier an- 
alytical results regarding settling in a quiescent disk 
l|Garaud fc Lin .2004) . However, the question of the effect 
of MHD turbulence was left unanswered. In this paper, we 
tackle this problem by means of local MHD simulations 
that use Eulerian finite difference codes. Our goal is to re- 
late the dust dynamics to the properties of the turbulence, 
and also to characterise the amount of settling as a func- 
tion of the particle size. Because of the complexity of this 



problem, we neglect dust coagulation and only consider 
populations of particles with a single specified size. 

The plan of the paper is as follows: in section El we 
review the characteristics of dust-gas interaction and de- 
scribe a simple model to account for the dispersion of a 
passive scalar in a turbulent medium with random veloc- 
ity field. In sectional we present the numerical methods 
we used together with the basic disk model which exhibits 
sustained MHD turbulence. We go on to study of the dif- 
fusion of very small dust particles in sectional Their evo- 
lution is analysed in the context of the simple diffusion 
model discussed in section |21 We then show in sectional 
that larger dust particles start to decouple from the gas 
and settle toward the midplane of the disk. In this case 
the evolution is found to be well described using an ad- 
vection diffusion equation incorporating settling together 
with the same diffusion coefficient that applied to the very 
small particles. In section we investigate the effect a 
dead zone is likely to have on these processes by setting 
up non-ideal MHD calculations that result in a signifi- 
cant midplane centered dead zone. However, gas motions 
are effective at preventing complete settling for 10cm sized 
particles in this case too. Finally, we discuss our results in 
section H 

2. Basic equations 

We adopt a two fluid model of the dust and gas circulat- 
ing in a protoplanetary accretion disk. The basic equations 
are those of MHD applied to the gas component together 
with those of hydrodynamics applied to the dust compo- 
nent. The model allows the two fluids to have different flow 
velocities and consequently exchange momentum through 
drag forces. 

The basic equations for the gas component are the con- 
tinuity, momentum conservation and induction equations. 
In a frame rotating with angular velocity fl = ilk, with 
k being the unit vector in the fixed direction of rotation, 
here called the vertical direction and fl being the magni- 
tude of the angular velocity, these take the form 

^ + V.(H = 0, (1) 

-:^ + (v-V)v + 2nxv = VP+ (VxB)xB 

ot p iirp 

- ElE^a:^ ^ (2) 

p riip 

— = Vx{vxB -7]VxB) . (3) 

Here, P is the gas pressure, p is the gas density, $ is sum 
of the gravitational potential, here due to a central mass 
and the centrifugal potential — 0^|rxkp/2 with r being 
the position vector, v is the gas velocity and B denotes the 
magnetic field. The resistivity of the gas , 77, is taken to be 
nonzero only in section |B| where we investigate the effect 
of non-ideal MHD. The last term in equation (O is the 
acceleration of the gas due to drag on the dust component. 
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The drag force acting on a single dust particle of mass irip 
is Fdrag (sce below) and pd is the dust density. 
The equations for the dust are 



— + {vd-'V)vd + 2flxvd 



drag 



(4) 
(5) 



with Vd being the velocity of the dust component. 



2.1. The drag force 

The dust interacts with the gas through drag. In this pa- 
per, we consider only dust particles that are small enough 
in siz e that they are in the Epstein regime ijWeidenschillind 
Il977tl . The drag force Fdrag acting on a single particle 
then takes the form 



■ drag 



rn. 

T. 



-{V ~ Vd) 



(6) 



This drag force is proportional to v — Vd which is the 
relative velocity between the gas and dust components, 
and Ts, the dust stopping time, defined by 



Psa 

PCs 



(7) 



where a is the dust particle radius, Ps is the dust particle 
density and Cj is the local speed of sound which, as we 
shall take an isothermal equation of state, is independent 
of height. 

When is taken to be the Keplerian angular velocity 
at some disk radius, the size of the particles there can 
be measured in terms of the dimensionless parameter flTs 
through the relation 



(8) 



where H = c .^ /^ is the disk scale height. Following 
Ijohansen et al.l l)2004j) . when we take Ps/p = 10^" and 
H = 10^^ cm, equation (jHJ gives 



100 (nTs) cm. 



(9) 



The above relation is also appropriate for the standard 
minimum mass solar nebula at 5.2AU. 

Using an a pproa ch based on Boltzmann averaging, 
ICaraud et alJ ((2004') argued that the evolution of the 
dust particle distribution can be accurately determined 
by modelling it as a pressureless fluid as long as the di- 
mensionless parameter flTs is less than unity. This is the 
case for the dust particles over most of the disk domains 
considered in this paper. Accordingly we adopt the two 
fluid description in preference to, for example a much more 
computationally demanding N-body approach. We also 
assume throughout this paper that the dust particles re- 
main electrically neutral so that they do not react to the 
magnetic field. 



2.2. Single fluid with advection diffusion equation and 
settling 

We here note that in the limit Ts 0, the two fluid de- 
scription can be reduced to one of a single fluid with MHD 
together with an advection diffusion equation describing 
the evolution of the dust. To see this we introduce a den- 
sity p and a velocity v associated with the combined fluid. 
These are defined through 



p = p + Pd 
and 

pV = pV + pdVd. 



(10) 



(11) 



Adding the continuity equations ^ and Q and eliminat- 
ing the drag force by combining equations Q and (jSJ and 
subsequently letting ^ while assuming v—Vd = 0{ts) 
gives 



dv 



0, 



-VP 



1 

47r 



(12) 

(VxB)xB 

(13) 



Correct to first order in Ts, we may set Vd — v in the 
left hand side of equation jSJ and then use that equation 
to find Vd located on the right hand side. Substituting this 
into equation Q and after making use of equation H13|l 
gives rise to the equation 



^ + V.(P.«)=0 
with 



where 



PTs_ 



^ -VP 



— {VxB)xB 



(14) 



(15) 



(16) 



Using the above in equation (|14|l and setting pd ~ Xjj, 
p = (1 — X)p, results in the following equation for the dust 
mass fraction X 



-p(^+v-Vx] + 



V • 



dt 



XpTs 



(l-X)V(pc2) , (Vxi?)xB 



P 



- V • {pclTsXVX), 



Airp 



(17) 



This can be interpreted as an advection diffusion equation 
where the advection velocity Ua is the sum of the mean 
fluid velocity and a settling velocity such that 



P 



(l-X)\/{pcl) ^ {VxB)xB ' 
P 47rp 



(18) 
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We remark that for small X and no magnetic field, the 
vertical component of the settling velocity is — with g 
being the vertical acceleration due to gravity. In addition 
there is a small diffusion coefficient arising through the 
second derivatives of X with respect to the coordinates 
such that D = pXcj.Ts. In practice, when turbulence is 
present, this is overwhelmed by the effect of the mean 
fluid velocity, that being particularly the case when X is 
small. 

2.3. A simple theory for dust diffusion 

The above analysis suggests that the dust mass fraction 
is advected as a passive scalar. When the gas is turbulent 
part of the velocity field will be turbulent. The charac- 
teristic correlation time associated with the turbulence is 
expected to be a fraction of an orbital period (see below), 
whereas the underlying settling takes place on a signifi- 
cantly longer timescale ~ l/(f2^rs). Thus we separate the 
effects of turbulent diffusion and settling and consider the 
simplest case of the dispersion of a passive scalar in steady 
state, homogeneous and isotropic turbulence. This means 
that effects related to the imperfect coupling between dust 
and gas such as turbul ent enhancement of gr ain-grain 
colhsions are negl ected (TVoelk et all 11980 : iMorfill. 1985.) . 
It is well known llTavloTlh 9211: lRa,tcheloT<ll 949h that the 
mean square displacement of the scalar from its location 
at some initial time can be expressed in terms of the fluid 
velocity fluctuations. We adopt a Lagrangian approach 
and for simplicity consider only the vertical z-direction. 
Extension to consider the other coordinate directions is 
straightforward but uninformative in our case. 
Let z(zcj, t) be the position of a particle which was located 
at z = zq at i = 0. Its displacement Z{zq, t) is given by 



Z{zo,t) ^ z{zo,t) - zq = / U{zo,t')dt' , 

Jo 



(19) 



where U{zQ,t) is the vertical component of the velocity 
of the particle as deflned in a Lagrangian sense as seen 
along its path. It relates to its Eulerian equivalent Vz{z,t) 
through U{zo, t) = Vz{z, t). We are interested in the mean 
square deviation of particles from their initial positions. 
Using the fact that U{zo,t) and Z{zo,t) are related by a 
time integral, the derivative of Z^(zo,t) with respect to t 
can be expressed in terms of the velocities as: 



dZ\zo,t) 
dt 



= 2 



U{zQ,t')U{zQ,t)dt' 



= 2 / Vz{zizo,t'),t')vz{z{zo,t),t)dt' . (20) 
Jo 

To obtain a mean square particle displacement we take 
an ensemble average over many particle realisations, indi- 
cated by angular brackets, to obtain 

d < Z2(zo,i) > 



dt 



2 / Szz{t,T)dT = 

Jo 



2/ <Vz{zizo,t-T),t-T)vz{zizo,t),t) > dr , (21) 
Jo 



where Szzit,T) is the velocity correlation function which 
depends only on the properties of the turbulence. For tur- 
bulence in a statistically steady state, as there is no pre- 
ferred time, it should depend only on the time difference 
t — t'~T. Thus any value of t may be adopted in Szz{t, r) 
so that 



Szz {t,T)^ Szz (i = 0, t) = Szz {t = t,t), 

which finally writes 

Szz{t,T) =< Vz{z{Z(3,t),t)Vz{zo,Q) > ■ 



(22) 



(23) 



When T — > 0, Szz becomes < >, i.e. it is a measure 
of the mean square turbulent velocity fluctuations. In the 
opposite limit , when t —^ oo, Szz tends to zero as the 
velocities become uncorrelated. 

The calculation of the velocity correlation described 
above is greatly simplified if the integral in H21I) can be 
performed at a fixed location using Eulerian velocity com- 
ponents, such that Vz{z{zq,t),t) is replaced by Vz{zo,t) 
with zq being fixed. This is possible when the distance 
moved by a particle during the turbulence correlation time 
is small compared to the c haracteristic scale of the turbu- 
lence (^Biferale et 811119951) as would always be the case if 
the turbulent velocities were reduced in magnitude while 
retaining their spatial and temporal characteristics. This 
situation seems to be a reasonable approximation in the 
case of the MRI where the turbulence is subsonic and we 
note that it would be even better in dead zones where 
the turbulent velocity fluctuations are reduced in magni- 
tude. In this case 5*22 simply scales down as the square 
of the velocity amplitude. Furthermore for the turbulence 
we consider, there is no preferred location or particles in 
the above ensemble averages and all fixed locations can 
be assumed to get complete information about all particle 
realisations. Then it is reasonable to replace the ensemble 
average over all possible particle realisations as a function 
of time, by one over all possible fixed locations using the 
velocities there at that time. Therefore, from now on, we 
adopt 



Szz{t) =< Vz{z, t)Vz{z, 0) > 



(24) 



where, the ensemble average is over fixed locations z and 
over different realisations obtained from a set of models. 
We also define the quantity Dt{t) to be 



Dt{t) 



Szz{T')dT'. 



(25) 



2.4. Diffusion coefficient 



As we shall see, for large r, the quantity D{t) defines an 
effective diffusion coefficient. For small r, Szz takes on a 
flnite value. Thus w e expect Dt( t) oc t. This is known as 



At 

the ballistic regime llTavloilll92l]) 

Dt{t) ^<v'l> t (for smaU r) 



(26) 



When r becomes large, 5*22 goes to zero, so we can postu- 
late that D(t) goes to a flnite limit, meaning that < Z^ > 
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would scale like t. This is the diffusive regime. If Tcorr is 
a typical correlation timescale for the turbulence, then we 
expect, from simple dimensional analysis 



Dt{t) 



(for large r). 



(27) 



However, at this point we insert a cautionary note. Strict 
convergence as r ^ oo would require ensemble averaging 
over an infinite number of realisations. In practical cases 
these will be finite in number and evaluating Dt{t) for 
increasing r corresponds to working out the average of a 
finite number of random walk displacements for increasing 
numbers of steps. Although this might appear reasonable 
for modest t, it would not be expected to show strict con- 
vergence for very large r in practical cases. For this reason 
a time span of 8 orbits equivalent to about 50 correlation 
times is adopted below. 

In this context we comment that the modelling of dust 
spreading using a diffusion coefficient obtained by the 
above procedure, on account of the averaging involved, of 
necessity only describes evolution occurring on time scales 
significantly longer than the correlation time associated 
with the turbulence. Behaviour occurring on the correla- 
tion time scale or faster cannot be meaningfully considered 
in this approach. 

3. Numerical Methods 

3.1. The algorithms 

In this paper, we use two well known finite^ differ- 
ence codes: ZEUS-3D fst one & Norman I9 92albl) and 
NIRVANA (|Zieg1er fc York Jl 9971) . Both are used here to 
solve the MHD equations in the shearing sheet approxi- 
mation (^Goldrcich & Lvndcn-Bell 1965), including verti- 
cal stratification. In this approximation a small region of 
the disk is considered in the neighbourhood of some point 
rotating in circular orbit with the local Keplerian angu- 
lar velocity. Local Cartesian coordinates are used with the 
X axis along the line connecting the origin to the central 
mass but pointing away from it, and the y axis directed in 
the direction of orbital motion. When dust is absent the 
equations solved are l^-El) with Fdrag = adapted to 
the shearing sheet by taking 



'Sn^xl 



(28) 



with i being the unit vector in the x direction. The 
evolution of the m agnetic field is calcula ted using the 
MOC-CT algorithm l|Hawlev fc Sto nc'l995V such that the 
solenoidality constraint on the magnetic field is satisfied 
at all times. 

In the two codes, we employ two different algorithms 
to describe the evolution of the dust particles that are the 
main focus of this paper. ZEUS-3D was extended to de- 
scribe the dust pa rticle evolution by m eans of a second, 
pressureless fluid l)Fromang fc Nelsonl[2005 ) using equa- 
tions H4I5I) together with appropriate drag forces in the 
low X limit. Because of the short stopping time of the 
dust particles, the effect of these is computed implicitly. 



On the other hand, NIRVANA solved the MHD equa- 
tions (|f 2I13() together with Consistent with the low X 
limit the mean flow velocity v was used in |j2Jl. The evo- 
lution of the dust mass fraction was calculated by solving 
equation H17|l but with the use of an advection velocity 
limiter. This was applied such as to restrict the settling 
velocity Ua — v to be less than O.lc^ in magnitude. This 
is required in the low density upper layers where dust gas 
coupling becomes weak and where the theoretical descrip- 
tion breaks down. However, because it applies only in the 
upper layers, calculations of settling dust are insensitive 
to this cut off that is required for numerical reasons. 

3.2. The model properties 

In this section, we describe the underlying disk model we 
used and highlight some of the properties of the MHD 
turbulence that are of importance for the dynamics of the 
dust that we go on to study. 

The initial disk set up in the ab s ence of dust is very 
similar to the model of IStone et a1.l (|l99fil) . The parame- 
ters are those of the standard shearing box. The equation 
of state for the gas is isothermal, P = pc^. Because of 
the vertical stratification, the initial density distribution 
is given by: 



p = p = poe 



(29) 



where pq is the midplane density and H — c^/fi is the disk 
scale height. The Cartesian box extends over the domains 
[-H/2,H/2] in X, [0, 27rif] in y and [-3H,3H] in z. We 
use the standard p eriodic boundary conditions in shearing 
coordinates in x llHawlev et all Il995t iBalbus fc Hawle\i 
EM) and periodic boundary conditions in y. One needs to 
be careful with the vertical boundary conditions. Because 
of the vertical density stratification, in particular the non 
vanishing vertical component of the gravitational force on 
the boundary, using simple periodic boundary conditions 
in that direction leads to unphysical fluctuations gener- 
ated in the box. 

In ZEUS-3D, we therefore made the vertical gravita- 
tional force continuous at the boundary by writing the 
gravitational potential as 



^^n'iz'-Hi),o 



(30) 



In practice, this means taking vertical gravity into account 
only in the region z e [—Hq,Ho]. We found that taking 
Ho = 2.AH produces satisfactory results. 

In NIRVANA, we further checked for the existence of 
numerically generated fluctuations by applying the addi- 
tional procedure of making the vertical gravitational ac- 
celeration continuous throughout. To do this the usual ac- 
celeration was applied for z S [— 2.25i/, 2.25i?]. This was 
then reduced to zero in ±[2.25i?, 2.4H] by means of linear 
interpolation and set to zero for \z\ > 2.AH. As indicated 
below, similar results were obtained with both codes. 
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time (orbits) 



Fig. 1. Time history of the total stress {dashed line), the 
sum of the Maxwell (upper solid line) and the Reynolds 
stresses (lower solid line) obtained with ZEUS-3D. All 
stresses are normalised by the midplane pressure. After 
the linear growth of the MRI, the total stress reaches a 
nonzero quasi-steady state, showing that the turbulence 
is sustained. 



Fig. 3. Gas distribution in the (r,z) plane. Density fluctu- 
ations across the disk are clearly visible, showing that the 
entire disk is turbulent at this stage. 




H4M 



time (orbits) 



Fig. 2. As in FigureHbut obtained with NIRVANA. 



In both cases the computational box is initially 
threaded by a magnetic field with zero net flux, being of 
the form 



B, = Bosin(2Trx/H) . 



(31) 



Bq is calculated such that /3, the ratio between the thermal 
pressure and the magnetic pressure is initially 100. Finally, 
the numerical resolution is (N^,Ny,N,_) = (32,100,192) 
for ZEUS and (36, 100, 196) for NIRVANA. Models with 
no dust are initiated by imposing a small random velocity 
component and subsequently run for ~ 100 orbits. 

We describe here simulation results, focusing on those 
obtained with ZEUS-3D. The results obtained using 
NIRVANA were very similar and accordingly only a few 
need to be illustrated here. As expected, the initial set up 




Fig. 4. Velocity dispersion as a function of height, nor- 
malised by the sound speed. The three curves corresponds 
to the radial (dashed line), azimuthal (solid line) and ver- 
tical (dotted line) velocity dispersions. 



is unstable to the MRI. Its early growth is illustrated in 
figure n Corresponding results obtained with NIRVANA 
are plotted in Figure |2] Both figures show the time his- 
tory of the volume averaged Maxwell stress (upper solid 
line) and Reynolds stress (lower solid line), normalised by 
the average midplane pressure Pq- They are respectively 
defined as 



rpMax 

nRey 



^dr. 

4tt 



rcf) 



p(Vx - Vx)(Vy - Vy)dT , 



(32) 
(33) 
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where Vx and Vy are respectively the verticahy and az- 
imuthally averaged radial and azimuthal velocities and dr 
denotes the element of volume. The dashed line in figure ^ 
is the sum of the Maxwell and Reynolds stresses: 

rpMax I rpRey 

a = '•^ I '"^ (34) 

After a few orbits, it peaks at a maximum as the flow 
breaks down into turbulence and attains a saturated state. 
For the remaining part of the simulation, a varies noisily 
between 0.006 and 0.03, consistent with previous studies 
ijStone et al.lll996(l . We remark that there are quite large 
fluctuations in the stresses. At the present time, their ori- 
gin is not clear. The origin of these large variations in our 
stratified zero net flux models should be addressed in the 
future using higher resolution calculations. In any case, 
the results show clear evidence that MHD turbulence is 
sustained for a very long time. The state of the disk is 
illustrated in flgure|31 where we plot a typical snapshot of 
the density in the (r,z) plane. Turbulent density fluctua- 
tions superposed on the overall vertical stratification are 
obvious from this figure. 

In section 12.31 we argued that the diffusion of small 
dust particles depends on velocity fluctuations in the disk. 
In figure ^ we plot their vertical profile normalised by 
the sound speed. The dashed curve corresponds to the ra- 
dial velocity fiuctuations, while the azimuthal and vertical 
fluctuations are respectively given by the solid and dotted 
curves. The midplane values given by this plot are 

(<5^;2)l/2^0.11c,, (<5«2)i/2 ^ q.OSc. , ^ o.07c, . 

These values are consis tent with numbers previously re- 
ported in the literature Ist one et aLlll996() . 

4. Vertical diffusion of small particles 

In this section, we present the results we obtained for very 
small particles using ZEUS-3D. For a given disk model, 
simulations of dust evolution are defined by the value of 
r^Tg. This is actually a function of gas density. To obtain 
a fixed parameter defining a given simulation, fir^ is eval- 
uated using the initial uniform midplane disk gas density. 
Where this is not stated explicitly it should be taken as 
read. We focus here on particles for which 17 Ts = 10^^ 
in the midplane of the initial disk. Equation @ indicates 
that this corresponds to micron size particles. This small 
value of JIts means that the dust particles are very well 
coupled to the gas and behave almost like a passive scalar. 

The dust is initially distributed in a thin layer around 
the disk midplane. Initially, the vertical profile of pd is 
taken to be a Gaussian with a thickness Hd = 0.2H. Under 
the effects of MHD turbulence, this small layer broadens 
with time. We found that its precise evolution depends 
largely on the particular time at which the dust particles 
are released in the disk. In order to improve the statis- 
tics, we computed nine simulations, by restarting the disk 
model described in section rOI at times t = 40, 45, 50, 55, 
60, 65, 70, 75 and 80 orbits. 



Fig. 5. The first four snapshots show the dust distribution 
in the (r,z) plane obtained with ZEUS-3D at times t — 3, 
5, 10 and 15 orbits after the dust is introduced (from left 
to right) . The last snapshot on the right represents the gas 
distribution at i = 15 orbits. At that time, the micron- 
sized dust particles are well mixed in the entire disk. The 
dust density distribution is also seen to be well correlated 
with the gas density. 



Figure illustrates the typical evolution of such a 
model. From left to right, the first four snapshots on the 
left show the dust density distribution in the (r,z) plane 
at time t = 3, 5, 10 and 15 orbits (measured after the dust 
has been introduced). The last snapshot shows the corre- 
sponding gas distribution at t = 15 orbits. The dust is seen 
to spread quite rapidly. At 15 orbits, it is almost filling the 
entire vertical extend of the disk. At later times, this dis- 
tribution does not change qualitatively. By comparing the 
dust and the gas distribution at t = 15 orbits, one can also 
see that they are very well correlated. We also found that 
the initial distribution of the dust has very little effect on 
this final state: starting with an initial dust distribution 
such that the gas-to-dust ratio is uniform gives an end 
product almost indistinguishable from the last two snap- 
shots of figure [S] 

In order to make a more quantitative comparison be- 
tween the results of these simulations and the simple the- 
ory presented in section 12.31 we computed the functions 
Szz and Dt according to equations H24() and (|25() . To cal- 
culate the former, we first volume averaged the velocity 
product in radius and azimuth. The vertical average is 
taken only within H around the disk equatorial plane. 
Indeed, the behaviour of the turbulence away from the 
midplane is likely to depart very much from being homo- 
geneous and isotropic and may be infiuenced by the ver- 
tical boundary conditions. To reduce the statistical noise, 
the result is then averaged over the nine different models 
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Fig. 7. Diffusion coefficient as a function of time, normalised by Csi?, obtained from the simulations done with ZEUS- 
3D. It is computed using equation H25|) and combining an average over a part of the computational box (|z| < H) and 
over a set of models (see the description in the text). As shown on the left panel, the diffusion coefficient is observed to 
saturate to a well defined value D ^ 5.5 x lO^^CgH after an initial rise. The dashed line plots the analytical estimate 
computed from equation (|27|l . The right panel shows an enlargement of the early evolution of Dt, compared with its 
expected early behaviour as given by equation l|26|) for small time (dashed line) 



0.025 



0.010 



time (orbits) 

Fig. 8. As in Figure [3 but for runs performed with NIRVANA. The lower curve applies for \z\ < H, and the upper 
curve for H < \z\ < 2H. In the former case the diffusion coefficient is observed to attain D ^ 5 x lO'^CsH after an 
initial rise in good agreement with the ZEUS results. In the latter case the estimated diffusion coefficient is less stable 
but is about 4 times larger. In that case the significance of the diffusion coefficient is less clear for the reasons given 
in the text. 



we computed. The resulting function Szz{t) is plotted in 
figure|n|using the solid line. The dotted line represents the 
time variation of the function 



Szzir) = (0.07c,)2e"*/^ 



(35) 



where Tcorr — 0.15 orbits. Given the good agreement be- 
tween the solid and dotted lines on figure El we conclude 
that Tcorr IS a measure on the typical timescale over which 
the fluid velocities become uncorrelated. 

Integrating Szz over time gives Dt{t). The left panel of 
Figure H shows the time evolution of Dt/{csH). Figure|Hl 
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0.0 0.5 1.0 1.5 2.0 
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Fig. 6. Time history of the function Szz{t), averaged be- 
tween the nine models that were run with ZEUS-3D. 
The solid Hne show the function as it is extracted from 
the models while the dotted line represents the function 
5*0 exp {—t/Tcorr), with 

Tcorr — 0.15 orbits. This dotted 
curve is seen to match nicely the numerical result. 



gives results of similar calculations done from runs with 
NIRVANA (see below) . Dt is first observed to rise at early 
times. This is in good agreement with equation H2t)|) . In 
fact, the prediction given by equation (|26|l is represented 
by the dashed line on the right panel of figure which 
is an enlargement of the left panel at small times. It uses 
iSvlY/"^ — O.OTc^, as obtained in section IT^ This early 
lin ear rise of the di ffusion coefficient was also observed 
bv ICarballido et all |2005). Here, we see that it is natu- 
rally understood in terms of the fluid velocity correlations. 
After this initial rise, Dt{t) is observed to reach a roughly 
constant value of 5.5 x lO^^CgH . This value nicely com- 
pares with the naive estimate of equation (|27|l . Indeed, 
taking the value of 5v1 and Tcorr derived above, one obtain 
Dt ^ 4.6 X IQ^^CgH , a value which is surprisingly close 
to the result of the numerical simulations. This analytical 
estimate is represented by the dashed line on figure [7| 

If the dust distribution undergoes diffusive evolution 
as indicated by the simple theory, pd should satisfy a dif- 
fusion equation of the form 



■ d 






m 



(36) 



where D is a constant diffusion coefficient. We tested this 
hypothesis by solving equation (|36l) . To do so, we took 
an initial dust distribution identical to that of the numer- 
ical model and investigate three different values for the 
diffusion coefficient D. First, we used D/{csH) = Dth = 
5.5 X 10"'^, i.e. the value given by our simple theoretical 
model in terms of the velocity fluctuations of the fluid. 



If this theory is correct, we expect this solution to be 
close to the numerical solution. In order to check its accu- 
racy, we also computed the solution of equation H36|l using 
D/{csH) = 10-3 and D/{c,H) = 10'^. 

The comparison between the analytical solution of 
equation H36|l and the results of the numerical simula- 
tions is shown in figure |H| It is illustrated by nine pan- 
els. They correspond, from top left to bottom right, to 
times t = 0.48, 0.80, 1.12, 1.44, 1.76, 2.08, 2.40, 2.72 and 
3.04 orbits (measured after the dust was introduced into 
the disk). On each panel, there are four curves. The solid 
line shows the vertical profile of the dust-to-gas ratio, av- 
eraged between the nine models we ran and normalised 
by the value in the equatorial plane at t — 0. The three 
dashed curves are the solutions of equation using the 
three diffusion coefficients mentioned above. Of course, the 
smaller the value of D, the less the dust is spread across 
the disk at a given time. 

The agreement between the numerical results and the 
simple model is fairly good. In all panels, the solid line is 
seen to have the best agreement with the middle dashed 
curve, calculated using D — Dt- The agreement is espe- 
cially good at low altitude (typically z < H). This was to 
be expected: first, this is the region in space where we per- 
formed the volume average used to calculate the function 
Szz- Second, the theory supposes that the turbulence is 
homogeneous and isotropic. This results in only possibly 
reasonable modelling near the midplane. The anisotropy 
of the turbulence is expected to increase away from the 
midplane, as the density stratification becomes stronger. 
From figure El there are some indications that the es- 
timated diffusion coefficient might increase with height. 
Indeed, at late times, for z > H, the solid line shows 
better agreement with the dashed curve computed using 
D/{csH) — 10^^ than at lower altitudes, indicating that 
dust particles are spread more efficiently than the simple 
theory suggests. This is also supported by results of runs 
performed with NIRVANA to test the simple diffusion the- 
ory. As with ZEUS-3D we considered runs restarted from 
the basic disk model after 33, 40, and 60 orbits. We evalu- 
ated the function Dt{t) by performing ensemble averages 
both over grid points such that \z\ < H, and for grid points 
such that H < \z\ < 2H and then averaging over the mod- 
els. Results are shown in Figure |S1 The indication is, in 
agreement with ZEUS, that for \z\ < H, Dt approaches 
- 5 X lO-^c^iJ, but for if < \z\ < 2H, Dt is 4. times 
larger. 

5. Vertical settling of larger particles 

In the previous section, we have studied the effect of MHD 
turbulence on very small dust particles, which behave es- 
sentially like passive scalars. In this section, we will fo- 
cus on larger particles for which the settling processes 
described in the introduction arc fast in the absence of 
turbulence. 

We focus on three sizes, for which the parameter 
is respectively 0.001, 0.01 and 0.1 in the midplane of the 
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Fig. 9. Vertical profile of the dust to gas ratio as a function of time. The different panels correspond, from top to 
bottom and from left to right, to t = 0.48, 0.80, 1.12, 1.44, 1.76, 2.08, 2.40, 2.72 and 3.04 orbits. On each panels, the 
solid line is the result of the simulations performed with ZEUS-3d, averaged over the nine models described in the 
text. Three dashed curves are plotted. They show the solution of equation H36() . computed using D/icgH) = lO^'^, 
D/{csH) — 5.5 X 10^^ and D/{csH) = lO^'^. The smaller the value of D, the smaller the diffusion. 



initial disk. Using equation (j^jl, these values correspond 
respectively to 1 mm, 1 cm and 10 cm. In the absence of 
turbulence, the typical settling timescale can be written 
I Dullemond fc Dominikll2004,l 



where Torb is the orbital time. In our three cases, U,Ts — 
0.001, 0.01 and 0.1 correspond to Tsett = 160, 16 and 1.6 
orbits respectively. Thus in the absence of turbulence, set- 
tling would be important within the duration of the sim- 
ulations in the last two cases. 

At t = 40 orbits, we introduced the particles into the 
turbulent disk model with the same vertical distribution 
as for the small particles discussed in section^J Their evo- 
lution was then followed until the vertical profile of the 
dust-to-gas ratio reaches a steady state. Once again, by 
running models in which the initial dust-to-gas ratio is 



uniform, we checked that this final distribution does not 
depend on the initial conditions. 

The results obtained with ZEUS-3D are illustrated in 
figure [TUl This shows the typical distribution of the dust 
particles in the (r,z) plane after they have reached a steady 
state. From left to right, the snapshots corresponds to 
flTs = 0.001, 0.01 and 0.1. While the smallest dust par- 
ticles are still efficiently spread across the disk (compare 
with the fourth panel of figure [SJ, there are some signs of 
vertical settling for the larger particles. This is obvious in 
particular for the largest particles. But note however that 
turbulence is quite efficient at preventing the complete 
settling of these particles. As discussed above, a very thin 
dust subdisk would form in just a few orbits in a quiescent 
disk. We have confirmed that very similar results for the 
degree of settling as a function of flTg (evaluated in the 
midplane) are obtained with NIRVANA which employed 
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Fig. 10. Dust distribution in the (r,z) plane, obtained 
with ZEUS~3D, after a quasi-steady state has been 
reached. From left to right, the different panels correspond 
to = 0.001, 0.01 and 0.1. 




Fig. 11. As in Figure ^1 but for results obtained with 
NIRVANA 



an advection diffusion treatment of the dust mass fraction. 
Steady state distributions corresponding to those shown 
in Figure ^| are plotted in Figure ^2 Again these are 
found to be initial condition independent for distributions 
initiated in the midplane regions. 

Once again, we can make use of the simple the- 
ory developed in section 12.31 to interpret these results. 
However, in this case the full advection diffusion equa- 
tion H17|) incorporating settling must be used with an 
anomalous diffusion coefficient. Neglecting Lorentz forces. 



0.1000 



0.0100 r 




Fig. 12. Steady state vertical profile of the dust-to-gas 
ratio when flxs — 0.01 {solid line). It has to be compared 
with the vertical profile calculated using equation (|36() . 
shown with the dotted line for three different values of the 
dimensionless diffusion coefficient D/{csH): 10~^, 5.5 x 
10-3 and 10~2. 



0.1000 




Fig. 13. Same as figure IT^ but for the case ^Its = 0.1. 



in the low X or Pd/p limit t his gives llDubrulle et alJll995l: 
ISchrapler fc HennindlioO^ . 

^-^(.f^r.p,) = i?-^p-(^-jJ . (38) 

Assuming a steady state, this equation gives a simple first 
order ordinary differential equation for the vertical dust 
mass fraction profile. We calculate solutions for the same 
three values of the diffusion coefficient that were used in 
section 0] and compare them with the numerical simula- 
tions. The results are shown in figure lT^ for the case where 
= 0.01 and in figure IT!^ for the case rir^ = 0.1. We 
do not consider the case flTg = 0.001 because the dust- 
to-gas ratio is almost uniform at the end of the simula- 
tion. In these figures the solid line shows the steady-state 
dust-to-gas ratio obtained at the end of the simulations 
while the three dashed curves are the corresponding solu- 
tions of equation calculated with D/{csH) — 10^^, 
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Fig. 14. Same as figure IH but fo r the "larger dead zone" 
model of iFleming fc Stond l)2003(l . 



5.5 X 10^^ and 10^^. As in section^] there is a good agree- 
ment between the solid line and the middle dashed curve, 
for which the magnitude of the diffusion coefficient corre- 
sponds to that estimated using the velocity fluctuations 
of the underlying model. The two other values of the dif- 
fusion coefficient can clearly be ruled out. There is also 
a marked difference between the solid and middle dashed 
curve in the upper layers of the disk. This was also ob- 
served in figure El As pointed out in section^ this is due 
to the increase of the diffusion coefficient at high disk al- 
titudes and also to the fact that our simple theory breaks 
down because the turbulence ceases to be something that 
can be profitably modelled as homogeneous and isotropic 
at these locations. 



6. The effect of a dead zone 

All the results presented so far in this paper suppose that 
gas and magnetic field are perfectly coupled throughout 
the entire vertical extent of the disk. However, protoplan- 
etary disks are probably cold and dense enough that this 
perfect coupling is unlikely, at least in some regions of the 
disk. This situation has ge nerated the "layered accretion" 
paradigm l)Gammielll996l) in which the gas is turbulent 
only in the upper layers of the disk, while a quiescent (or 
"dearf") zone exists around the midplane. The extent of 
this dead zone is very uncertain, since it depends on the 
ionising source and on the chemistry | ^no et al. 200^ 
iFromang et"ai]l2002l: lllgner fc Nelsonll2005|) . but its exis- 
tence is likely. In this section, we investigate how the pres- 
ence of a dead zone would infiuence the results described 
in the previous sections. 

The problem of layered turbulence in local numerical 
simulations of stratified disks has already been studied by 
iFleming k, Ston3 l)2003j) . In this section, we reproduce one 
of their disk models by allowing the resistivity 77 to be a 
function of position. We choose the vertical profile of rj 



Fig. 15. Steady state vertical profile of the dust to gas 
ratio in the ideal MHD case {dashed line) and when a 
dead zone is present around the equatorial plane of the 
disk {solid line). The parameter fir^ equals 0.01 in that 
case. 



such tha t our model is identical t o the "larger dead zone" 
model of IFleming fc Stond l|2003|) : 



77 = 770 exp 



(39) 



where Yiq/Y,cr = 30 and ryo is chosen such that the 
Reynolds magnetic number Rem, defined by 



77 



(40) 



equals 100 in the midplane of the disk. This disk model was 
run for 100 orbits. The evolution is the same as that found 
by [Fleming fc Stone (2003 ). In particular, we found that 
density waves are excited in the dead zone by the turbu- 
lent motions in the active layers. The velocity dispersions 
of the three velocity components are shown in figure 1141 
as a function of z. This figure should be compared with 
figure^ Taken as a whole, there is a decrease of the ve- 
locity fiuctuations compared to the fully turbulent model. 
The presence of the dead zone is apparent through a de- 
crease in the azimuthal and vertical velocity fiuctuations. 
The latter have a root mean square dispersion of roughly 
0.03cs, thereby showing a decrease by about a factor of 
two compared to the fully turbulent model. However, since 
it does not vanish, a nonzero diffusion coefficient can be 
expected. 

Using this new underlying disk model and introduc- 
ing dust particles, we used ZEUS-3D to recalculate the 
models described in section |S1 for which Q,Tg — 0.01 and 
Q.Ts = 0.1. After a few orbits, the dust distribution reaches 
a new equilibrium state. It is represented in figure ITsl ffor 
Q,Ts = 0.01) and figure^] (for Q.Ts = 0.1). In both cases, 
the solid line shows the vertical profile of the dust-to- 
gas ratio in the nonideal case with a dead zone, while the 
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are seen to match very accurately the sohd and dashed 
curves that are their numerical analogues. 

The same procedure was followed in the model for 
which Qts = 0.01. In that case, we found Hd = 0.77 
and Hd = 0.37. The dust-to-gas ratio vertical profiles 
derived using these two values are plotted in figure 1151 
with dotted lines. The agreement with the solid line is 
quite good. However, there is a poor agreement with the 
dashed line. This is because the dust is spread to higher 
altitudes in that case. The hypothesis that Ts is a constant 
which we used to derive equation H41|) and H42|) starts to 
break down, which explain the discrepancy with the nu- 
merical result. 



Fig. 16. Same as figureEl but for the case flTg = 0.1. 



dashed line corresponds to the fully turbulent case. As ex- 
pected, the thickness of the dust sub-disk is smaller in the 
former case. 

It is possible to compare these results to those ob- 
tained using the simple model presented in section lT^ To 
do so, we first note that close to the midplane, the stop- 
ping time Tg is nearly constant. A steady-state solution to 
equation (|38|l can be written in that case as 



Pd 
P 



where the dust scale height Hd is given by 

Hd = 



D 



(41) 



(42) 



We first focus on the case VtTg = 0.1. When the disk is 
completely turbulent (or "active"), we found in section^ 
that the value of the diffusion coefhcient was D/{csH) = 
5.5 X 10~^. Using equation H42|) . this would give a scale 
height 



~ active 

Hd = 0.23 



(43) 



for the dust sub-disk. Next, we seek an estimate of the 
dust sub-disk scale height in the presence of a dead zone. 
By combining equation H27() and equation H42() , Hd can be 
related to the velocity fluctuations by 



Hd oc (sviy/^ . 



(44) 



Given the smaller value we obtained for {Sv'^Y^^ in the 
"larger dead zone" model (see figure I14|l , we therefore ex- 
pect 



~ dead ~ active . 

Hd ^ Hd /2 ' 



0.11 



(45) 



The two dust distributions computed with equation (|41|l 

~ active - dead , , 

using Hd and Hd (corresponding respectively to 
the fully turbulent case and to the "larger dead zone" 
model) are plotted in figure ITCl using dotted lines. Both 



7. Discussion 

In this paper, we studied the effects of MHD turbulence 
on dust settling by means on local numerical simulations 
performed using ZEUS-3D and NIRVANA, being Eulerian 
MHD codes using finite differences. 

We first investigated the case of very small particles 
which are strongly coupled to the gas. Turbulent velocity 
fluctuations were found to cause an initially thin dust sub- 
disk to spread. The time evolution of the vertical profile 
for the dust-to-gas ratio can be well modelled by a dif- 
fusion equation, with a diffusion coefficient D that can 
be expressed in terms of turbulent velocity correlations. 
We found that a simple analytical estimate of D can be 
obtained in terms of the mean square amplitude of the 
velocity fluctuations 5v1 and their correlation time Tcom 
both of which are properties of the turbulence alone: 



D ^ {6vlf^, 



(46) 



Similarly to ^ Turner et al. I |2006^, we found an increase 
of the diffusion coefficient with disk height. We also no- 
ticed that the ensemble averages used to calculate it show 
weaker convergence when the upper layers of the disk are 
included. While the diffusive description of dust spreading 
seems to work well in the neighbourhood of the midplane, 
it is less accurate at disk heights exceeding a few scale 
heights. 

A standard procedure in this type of analysis is to de- 
termine the value of the Schmidt number Sc, defined as 
the ratio between the anomalous viscosity and the dif- 
fusion coefficient. The standard ap proach in dust diffu- 
sion modelling is to take 5*^ = 1 IS chrapler fc Hen nind 
l2004t lllgner et al.li2004: Dullemond fc Dominik 20041 . In 
non ze ro net flux local simula tions of radial dust dif- 
fusion, J^arballido ct al.' ('2005') found Sc = 11, while 
I.Tohanse n 1 ^ Klahr (2005) reported Sr = 1.5 for yertica. 1 
diffusion in zero net flux simulations. iTurner et al.l l(200(Tf) 
also reported a near unity Schmidt number in their calcu- 
lations. It is worth comparing these values to the Schmidt 
number we can derive from our simulations. By averag- 
ing the total stress represented in figure ^ between 20 and 
100 orbits, one obtains a = 1.54 x 10^^. This, together 
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with the value of the diffusion coefRcient obtained from 
the velocity fluctuations gives 



an intermedia te valu e bet ween the meas u res o f 
IJohansen et"an l)2005(l and ICarballido et all l)2005|) . 

However, it is important to stress here that the origin 
of a non zero diffusion coefficient is on account of the 
velocity fluctuations and not in the transport properties 
of angular momentum. 

When dust particles grow to centimeter sizes, we found 
that they start to decouple from the turbulence and settle 
towards the midplane. The steady state profile of the dust- 
to-gas ratio is well approximated by the solution of an 
advection-diffusion equation. Even for particles as large 
as 10 cm, we found that the dust sub-disk is significantly 
spread since its semi-thickness Hd equals 0.23H, while the 
settling timescale in a quiescent disk is very short in that 
case (1.6 orbits). We note however that radial migration 
is important for particles of this size and the interplay be- 
tween that migration and MHD turbulence in a stratified 
disk could lead to complex phe nomena, such as local en- 
hancement of the dust density ijFromang fc Nelsonllioosfl 
that might affect this picture. 

Because they are cold and dense, protoplanetary disks 
are unlikely to have adequate ionisation to be turbulent 
everywhere. Therefore we also investigated the effect of the 
presence of a dead zone around the midplane. As expected, 
we found thinner dust sub-disks in that c ase. Similarly to 
previous studies l)Fleming fc Stondl2003|) . we found the 
dead zone is able to maintain significant activity (excited 
by the turbulent velocity fluctuations of the active zone). 
This activity is able to prevent the complete settling of 
10 cm size particles. However, we want to emphasise that 
for computational reasons, our analysis was limited to a 
case in which the mass of the dead zone roughly equals 
the mass of the active zone. We expect our result to be 
modified in cases where the mass of the dead zone is much 
larger than that of the active zone and therefore only apply 
them to dead zones that involve a modest fraction of the 
local surface density. 

Nonetheless for conditions appropriate to a minimum 
mass solar nebula, the work presented here that consid- 
ered MHD turbulence, t aken together with that of eg. 
iGomez fc Ostrikeil l(2005l) indicates that gravitational in- 
stability of the dust layer is unlikely and that the for- 
mation of objects of planetesimal size may depend on 
phenomena such as densification in vortices operating to- 
gether with vertical settling. 

For practical reasons, we neglected grain growth in this 
work. This is an important simplification, as dust parti- 
cles are likely to grow at the same ti me as they settle 
toward the equat orial plane of the disk l|Cuzzi et al ]ll996; 
iDullemond fc Dominik 2005) . Simulations of the evolution 
of an entire dust population through grain growth and tur- 
bulent stirring, that take account of both radial and ver- 
tical disk structure, are very challenging with present day 



computational capabilities, but will have to be performed 
in the future. 
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